clear all
global Ux Uy  h  g  p  q r tmax Lx Ly L dx dy Nx Ny y C E F A;
%----- Physical Parameters----------

Ux = 1; % Current Velocity in the X direction. Assumed Constant
Uy=  1;  % Current velocity in the Y direction
h = 50; % height of the bed, assumed constant
g=  9.81; % Gravitational Constant


%------------ Simulation Parameters % -----------------
Lx=  100;         % Length in the X direction
Ly =100;            % Length in the Y direction
Nx = 8;         % number of axial nodes in the X direction
Ny=8;          % Number of Nodes in the Y direction

steps_x = Nx-1;    % number of axial steps
steps_y = Ny-1;
dx = Lx/steps_x;   % axial step spacing, delta x
dy = Ly/steps_y;
%% Elliptic Equation

cc =zeros(Nx+3,Ny+3);
ww =cc;
ss = cc;
nn=cc;
ee =cc;

%% The Inner Points
for ii=2:Nx+1
    for jj= 2:Ny+1
%        cc(jj,ii) = -2*h*(1/(dx*dx)+1/(dy*dy))+h/3;
        cc(jj,ii) = 0.01;
        ww(jj,ii) = -h/(dx*dx);
        ss(jj,ii) = -h/ (dy*dy);
        ee(jj,ii) = h/(dx*dx);
        nn(jj,ii) = h/ (dy*dy);
        
    end
end
%% North and South Boundary Points

for ii = 2:Nx+1
    %   coeff = nn(ii,Ny);
    cc(Ny+1,ii) = cc(Ny+1,ii) + nn (Ny+1,ii);
    nn (Ny+1,ii) =0;
    cc (2,ii) = cc(2,ii) + ss(2,ii);
    ss(2,ii) =0;
end

%% East and West Boundary Points
 
for jj = 2:Ny+1
    %   coeff = nn(ii,Ny);
    cc(jj,Nx+1) = cc(jj,Nx+1) + ee (jj,Nx+1);
    ee (jj,Nx+1) =0;
    cc(jj,2) = cc(jj,2) + ww (jj,2);
    ww(jj,2) =0; 
end

%generatespy_index(Nx,Ny, cc,ss,ww,ee,nn);
%pause

[ ccc,css,cww,cee,cnn,m_cne,m_cnw,m_cse, m_csw  ] = rrb_index(Nx,Ny,cc , ss, ww);
generatespy9_index(Nx,Ny, ccc,css,cww,cee,cnn,m_cne,m_cnw,m_cse, m_csw );

[L1,U]=lu(A);

isequal(A,A')


%index = (j-1)*Nx + i;

%ww = index-1. ee = index+1 .. ss = index-Nx  .. nn = index + Nx